MODULE sefmm_quadrature

   USE fmm_global_paras
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: sefmm_init_quadrature,        &
             get_Coulomb_quadrature,       &
             get_roots_and_weights

   REAL(REALK), TARGET, SAVE :: exponents(1000)
   REAL(REALK), TARGET, SAVE :: weights(1000)

   INTEGER(INTK), SAVE :: GL_nquad = 200
   INTEGER(INTK), SAVE :: GL_per_cell = 10
   REAL(REALK),   SAVE :: GL_tmax = 10000d0
   LOGICAL,       SAVE :: SE_harri_quad = .FALSE.
!   INTEGER(INTK), SAVE :: GL_nquad = 40
!   INTEGER(INTK), SAVE :: GL_per_cell = 40
!   REAL(REALK),   SAVE :: GL_tmax = 30.0d0

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE sefmm_init_quadrature(scheme)

      IMPLICIT NONE
      TYPE(scheme_paras), INTENT(IN) :: scheme

      GL_nquad = scheme%GL_nquad
      GL_tmax = scheme%GL_tmax
      GL_per_cell = scheme%GL_percell
      SE_harri_quad = scheme%SE_harri_quad

   END SUBROUTINE sefmm_init_quadrature

!-------------------------------------------------------------------------------

   SUBROUTINE get_Coulomb_quadrature(mode,rts,wts)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN) :: mode
      REAL(REALK), POINTER :: rts(:), wts(:)

      IF ( (mode == 0) .AND. SE_harri_quad ) THEN
         CALL get_Harrison_Coulomb_quadrature(rts,wts)
      ELSE
         CALL get_GL_Coulomb_quadrature(rts,wts)
      END IF

   END SUBROUTINE get_Coulomb_quadrature

!-------------------------------------------------------------------------------

   SUBROUTINE get_GL_Coulomb_quadrature(rts,wts)

      IMPLICIT NONE
      REAL(REALK), POINTER :: rts(:), wts(:)
      INTEGER(INTK) :: i

      rts => exponents(1:GL_nquad)
      wts => weights(1:GL_nquad)

      CALL get_roots_and_weights(GL_nquad,GL_per_cell,zero,GL_tmax,rts,wts)

      DO i = 1, GL_nquad
         wts(i) = wts(i) * two/ROOTPI
      END DO

   END SUBROUTINE get_GL_Coulomb_quadrature

!-------------------------------------------------------------------------------
! Optimized Gaussian expansion of the Coulomb operator with precision
! of at least 10e-8 on [10e-5,1]

   SUBROUTINE get_Harrison_Coulomb_quadrature(rts,wts)

      IMPLICIT NONE
      REAL(REALK), POINTER :: rts(:), wts(:)

      INTEGER(INTK), PARAMETER :: nquad = 52
      INTEGER(INTK) :: i

      exponents(1 ) = 1.3372524762366868d11 
      exponents(2 ) = 8.0467192739036288d10 
      exponents(3 ) = 4.8419958241514414d10 
      exponents(4 ) = 2.9136009117055956d10 
      exponents(5 ) = 1.7532174747530922d10 
      exponents(6 ) = 1.0549735552210995d10 
      exponents(7 ) = 6.3481532107900659d9  
      exponents(8 ) = 3.8199111373359423d9  
      exponents(9 ) = 2.2985774753310111d9  
      exponents(10) = 1.3831365359548369d9  
      exponents(11) = 8.3228290858002536d8  
      exponents(12) = 5.0081451937458747d8  
      exponents(13) = 3.0135809077331903d8  
      exponents(14) = 1.8133799421750353d8  
      exponents(15) = 1.0911758996108682d8  
      exponents(16) = 6.5659977171864032d7  
      exponents(17) = 3.9509969363849716d7  
      exponents(18) = 2.3774569471066599d7  
      exponents(19) = 1.4306013528591281d7  
      exponents(20) = 8.6084429031350670d6  
      exponents(21) = 5.1800097407538342d6  
      exponents(22) = 3.1169981930574660d6  
      exponents(23) = 1.8756099387002403d6  
      exponents(24) = 1.1286219718397956d6  
      exponents(25) = 6.7913244132607723d5  
      exponents(26) = 4.0865842127987797d5  
      exponents(27) = 2.4590447345066979d5  
      exponents(28) = 1.4796956808832101d5  
      exponents(29) = 8.9038612357311147d4  
      exponents(30) = 5.3577736255235890d4  
      exponents(31) = 3.2239651392691467d4  
      exponents(32) = 1.9399758085236279d4  
      exponents(33) = 1.1673532360305863d4  
      exponents(34) = 7.0243843857770776d3  
      exponents(35) = 4.2268247930768600d3  
      exponents(36) = 2.5434325417535429d3  
      exponents(37) = 1.5304748626912268d3  
      exponents(38) = 9.2094178516074948d2  
      exponents(39) = 5.5416380390629165d2  
      exponents(40) = 3.3346029840785694d2  
      exponents(41) = 2.0065505753350410d2  
      exponents(42) = 1.2074136691414728d2  
      exponents(43) = 7.2654424320032992d1  
      exponents(44) = 4.3718781041547102d1  
      exponents(45) = 2.6307163144706104d1  
      exponents(46) = 1.5829948635381633d1  
      exponents(47) = 9.5249315234124400d0  
      exponents(48) = 5.7220041706777604d0  
      exponents(49) = 3.3624223407094093d0  
      exponents(50) = 1.7537139460449947d0  
      exponents(51) = 6.4705932650658966d-1 
      exponents(52) = 7.2765905943708247d-2 
      
      weights(1 ) = 1.0479544277680031d5
      weights(2 ) = 8.1291445822243042d4
      weights(3 ) = 6.3059049364932868d4
      weights(4 ) = 4.8915904045532969d4
      weights(5 ) = 3.7944848401804876d4
      weights(6 ) = 2.9434429047325397d4
      weights(7 ) = 2.2832762205449004d4
      weights(8 ) = 1.7711743950151233d4
      weights(9 ) = 1.3739287867104177d4
      weights(10) = 1.0657789571075259d4
      weights(11) = 8.2674214105396183d3
      weights(12) = 6.4131739752013645d3
      weights(13) = 4.9748040283865428d3
      weights(14) = 3.8590369818855305d3
      weights(15) = 2.9935182449329915d3
      weights(16) = 2.3221211966811754d3
      weights(17) = 1.8013075096471964d3
      weights(18) = 1.3973037965981704d3
      weights(19) = 1.0839114914325070d3
      weights(20) = 8.4080793916920919d2
      weights(21) = 6.5222852436674942d2
      weights(22) = 5.0594437698350613d2
      weights(23) = 3.9246936231794106d2
      weights(24) = 3.0444493025732431d2
      weights(25) = 2.3616293284245360d2
      weights(26) = 1.8319546607860352d2
      weights(27) = 1.4210773218655147d2
      weights(28) = 1.1023530215723992d2
      weights(29) = 8.5511334670538226d1
      weights(30) = 6.6332546980669662d1
      weights(31) = 5.1455246335384137d1
      weights(32) = 3.9914679842944927d1
      weights(33) = 3.0962472840916210d1
      weights(34) = 2.4018098812215013d1
      weights(35) = 1.8631233802429659d1
      weights(36) = 1.4452554123315050d1
      weights(37) = 1.1211083651910594d1
      weights(38) = 8.6966217584849718d0
      weights(39) = 6.7461123616573196d0
      weights(40) = 5.2330701805752999d0
      weights(41) = 4.0593785050153956d0
      weights(42) = 3.1489265907663571d0
      weights(43) = 2.4426740821107160d0
      weights(44) = 1.8948224052285526d0
      weights(45) = 1.4698450590705008d0
      weights(46) = 1.1401926133052701d0
      weights(47) = 8.8479121742292529d-1
      weights(48) = 6.9268638708061648d-1
      weights(49) = 5.8524457689702328d-1
      weights(50) = 5.7618252254532759d-1
      weights(51) = 5.9668881738899718d-1
      weights(52) = 6.0787990115110877d-1       

      rts => exponents(1:nquad)
      wts => weights(1:nquad)

      DO i = 1, nquad
         rts(i) = SQRT(rts(i))
      END DO

   END SUBROUTINE get_Harrison_Coulomb_quadrature

!-------------------------------------------------------------------------------
! First divide interval into logarithmic cells and then apply normal
! GL quadrature within cells

   SUBROUTINE get_roots_and_weights(npts,per_cell,x1,x2,rts,wts)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: npts, per_cell
      REAL(REALK),   INTENT(IN)  :: x1, x2
      REAL(REALK),   INTENT(OUT) :: rts(:), wts(:)

      REAL(REALK)   :: lo, hi
      INTEGER(INTK) :: n, nlo, nhi, ncell

      ncell = npts/per_cell
      lo = x1
      hi = x2/(two**(ncell-1))
      DO n = 1, ncell
         nlo = (n-1)*per_cell +1
         nhi = MIN(n*per_cell,npts)
         CALL get_GL_rts_and_wts(per_cell,lo,hi,rts(nlo:nhi),wts(nlo:nhi))
         lo = hi
         hi = two*hi
      END DO

   END SUBROUTINE get_roots_and_weights

!-------------------------------------------------------------------------------
! Return roots and weights of the Gauss-Legendre n-point quadrature 
! formula, given lower,x1 and upper,x2 integration limits

   SUBROUTINE get_GL_rts_and_wts(n,x1,x2,roots,weights)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: n
      REAL(REALK),   INTENT(IN)  :: x1, x2
      REAL(REALK),   INTENT(OUT) :: roots(:), weights(:)

      REAL(REALK)   :: xm,xl, z,z1, pp,p3,p2,p1
      INTEGER(INTK) :: m,i,j

      m = (n+1)/2
      xm = 0.5d0*(x2+x1)
      xl = 0.5d0*(x2-x1)

      DO i = 1, m
         z = COS(PI*(i-0.25d0)/(n+0.5d0))

         DO ! until z converged
            p1 = 1d0
            p2 = 0d0
            DO j = 1, n
               p3 = p2
               p2 = p1
               p1 = ((2d0*j-1d0)*z*p2 - (j-1d0)*p3) / j
            END DO
            pp = n*(z*p1-p2) / (z*z-1d0)
            z1 = z
            z = z1 - p1/pp
            IF ( ABS(z-z1) < 1d-15) EXIT
         END DO

         roots(i) = xm - xl*z
         roots(n+1-i) = xm + xl*z
         weights(i) = 2d0*xl / ((1d0-z*z)*pp*pp)
         weights(n+1-i) = weights(i)

      END DO

   END SUBROUTINE get_GL_rts_and_wts

!-------------------------------------------------------------------------------

END MODULE sefmm_quadrature
